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Abstract 



We introduce an operator description for a stochastic sandpile model 
with a conserved particle density, and develop a path-integral repre- 
sentation for its evolution. The resulting (exact) expression for the 
effective action highlights certain interesting features of the model, for 
example, that it is nominally massless, and that the dynamics is via 
cooperative diffusion. Using the path-integral formalism, we construct 
a diagrammatic perturbation theory, yielding a series expansion for the 
activity density in powers of the time. 
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I. INTRODUCTION 



Sandpile models were introduced some fifteen years ago as examples of self- 
organized criticality (SOC), or scale- invariance in the apparent absence of control 
parameters Subsequently the appearance of such "spontaneous" criticality 

was shown to result from a control mechanism that forces the system to a critical 
point marking a phase transition to an absorbing state @,^]. In fact, sandpiles with 
the same local dynamics as the original "self-organized" versions, but with strictly 
conserved particle density, p, exhibit an absorbing-state phase transition as p is 
varied. Thus the particle density is the temperature-like control parameter for these 
models. 



Sandpiles with strictly conserved particle density (so-called fixed-energy sandpiles 
or FES ||), exhibit absorbing-state phase transitions P-pH, which have attracted 
much interest of late, in connection with epidemics catalysis [|II||14[], and the 
transition to turbulence [|T5| — . While continuous phase transitions to an absorbing 
state fall generically in the universality class of directed percolation (DP) fi8| , |T9| , 
numerical evidence indicates that this is not so for sandpiles p0|-p3|. The non-DP 



nature of the transition has been attributed to the coupling of the order parameter 
(the density of active sites), to a second, conserved field (the local particle den- 
sity), which relaxes diffusively in the presence of activity, and remains frozen in 
its absence 
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Independent of the connection with SOC, fixed-energy sandpiles 
furnish intriguing examples of phase transitions far from equilibrium, incorporating 
cooperative relaxation in the form of activated diffusion. Until now, all quantita- 
tive results for FES have been numerical, based either on simulations |]20| , ^ , [2ll . |22| . 
or on coherent- anomaly analysis of a series of n-site cluster approximations |f23 . 
Phenomenological field theories for sandpiles have been proposed fl25|-f27|1 , but their 
analysis is far from straightforward. It is therefore of great interest to develop the- 
oretical approaches for FES. 



This paper is the first of a series analyzing a stochastic sandpile using operator 
and path-integral methods. In this work we establish the path-integral represen- 
tation, and use it to develop a time-dependent perturbation theory that yields an 
expansion for the activity density in powers of time. In subsequent work the series 
will be extended and analyzed. 



Our analysis uses two principal tools. The first is an operator formalism for 



Markov processes, of the kind developed by Doi [25|, and applied to various mod- 
els exhibiting nonequilibrium phase transitions |29| p3|. The second is an exact 
mapping, due to Peliti, from the Markov process to a path-integral representation 
|34] , |35|1 . This approach is often used to generate an effective action, which may then 
be analyzed using renormalization group (RG) techniques [j36| -|38ll- In the present 



instance, however, we use the formalism not as the basis for an RG analysis, but to 
generate a series expansion for the order parameter. The expansion variable is the 
time; the coefficients are polynomials in the particle density, p. 
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The balance of the paper is organized as follows. In Sec. II we define the 
model and the master equation governing its dynamics using an operator formalism, 
which is then mapped to the path-integral representation. Sec. Ill develops the 
perturbative expansion of the activity density in powers of time. In Sec. IV, we 
study some general properties of the diagrammatic expansion, allowing us to simplify 
and extend the analysis. A brief comparison with simulation results is presented in 
Sec. V. A summary and discussion of our results is given in Sec. VI. 



II. EVOLUTION OPERATOR 



version 



We consider Manna's stochastic sandpile in its fixed-energy (particle-conserving) 



The configuration is specified by the occupation number at 
each site; sites with n > 2 are said to be active, and have a positive rate of toppling. 
When a site topples, it loses exactly two particles ("grains of sand"), which move 
randomly and independently to nearest-neighbor (NN) sites. In this work, we adopt 
a toppling rate of n(n— 1) for a site having n particles. This choice of rate represents 
a slight departure from the examples studied previously, in which all active sites 
have the same toppling rate. The present rate leads to a much simpler evolution 
operator, and, on the other hand, should yield the same scaling properties, since 
sandpiles, like critical phenomena in general, exhibit a high degree of universality. 
(Close to the critical point, the density of sites with n > 3 particles is quite low, 
so that in practical terms our choice of rate should not alter quantitative properties 
greatly. Since studies of restricted-height sandpiles p2] , |2^| reveal that they belong to 
the same universality class as their unrestricted counterparts, there is good reason 
to expect that a small change in transition rates, that does not modify the sym- 
metry or conservation laws of the model, will have no effect on critical exponents.) 
Preliminary simulation results [[|l]] indicate that the model studied here exhibits 
a continuous phase transition at p c = 0.9493, in one dimension; the corresponding 
value for the model with the same toppling rate for all active sites is p c = 0.94885 



For simplicity we begin the analysis in one dimension; the generalization to d 
dimensions is straightforward. We represent the dynamics of this continuous-time 



Markov process via the master equation, written in the form p4"|,|3"5 



dt 



L|¥), 



where 



l*> = £KKM)IM> 



{rit} 



is the probability distribution. Here p({rii}, t) is the probability of the configuration 
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with occupation numbers {n«}, and the state is a direct product of states 

\n,j), representing exactly nj particles at site j. 

The evolution operator L is written in terms of creation and annihilation oper- 
ators, defined via the relations: 

a^rii) = riilrii-l) 

and 

Tii\rii) = |n» + l>. 

The evolution operator for the one- dimensional stochastic sandpile takes the form, 

1 



'I i \2 2 

-(7Ti_i +7Tj + iJ -TTi 



(1) 



It is readily seen that L conserves the number of particles. 

It is convenient to work in the Fourier representation; on a lattice of N sites with 
periodic boundaries, we introduce the discrete Fourier transform: 

a fc = , (2) 



with inverse 

1 



a ; 



1 N 



E e ^' ( 3 ) 



k 



(and similarly for 7tk, etc.), where the allowed vales for the wave vector are: 

2tt 2tt 2tt 2tt 

k = — 7T, — 7TH ,0, , 7T . (4) 

JV' N N N K J 

(To avoid heavy notation, we indicate the Fourier transform by the subscript k; the 
subscript j denotes the corresponding variable on the lattice.) In Fourier represen- 
tation, 



L = -N 3 E cj felifc2 7rfe 1 7rjt 2 afc 3 a_fe 1 _fe 2 _fc3 , (5) 



where LUk u k 2 — 1 — cosfci cos &2- 

Since L is in normal form (all operators it to the left of all operators a) , we can 
immediately write the evolution kernel, following Peliti's prescription |]34],|35|: 
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(6) 



fci,te,fc3 ) k 



with the boundary conditions ip k (0) = ( k and 4> k (t) = z k . (The dot denotes a time 
derivative; {z} stands for the set of parameters z k associated with each wave vector, 
and similarly for {(}■) The functional integrals are over the variables ip k and ip k : 

[v$Vi> = l[[v4> k [vi; k (7) 

k 



[The product is over the first Brillouin zone, Eq.(|]).] In the <i-dimensional case, Eq. 
(§) remains valid if we let 

Wk lf k a = 1 - A d (ki)A d (ki), (8) 

with 

1 d 

Ad(k) = - ^ cosfc Q . (9) 



a=l 



The wave vectors now range over the first Brillouin zone in d dimensions. 

The kernel is used to evolve the probability generating function, Q t ({zj}) = 
T,{ nj }P{{ n 3}^)T[jzf : or its Fourier equivalent, <$>t{{zk}) = J2{n k }P({nk},t) U k z k k , 

via a, 

**({**}) = 11/ ^^e-^^U t ({z k }, {Ck})M{<}) (t > 0). (10) 



We consider an initial product-Poisson distribution, with each site having the same 
mean number of particles, p, so that, 

$ ({^}) = exp[p^-l)], (11) 

3 

corresponding to 

%({z k }) = exp[Np(z k=0 -l)] , (12) 

in /c-space. Noting that 
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2tt 



/(C)e^^O = /(p), (13) 



we see that for a Poisson initial distribution, 

* t ({*}) = e-^l7t(M;C* = iVp**,o) • (14) 

Thus the generating function at time t is e~ Np times the r.h.s. of Eq. (|), evaluated 
with ^ fc (0) = Np5 kfi . 

The simplest observables of interest are the mean particle number (at site j), 

d§t({Zn}) 



(n,-) 



(15) 

z„=l 



and the mean activity, 

(16) 



d^ 2 



Zn=l 



In fc-space the condition Zj = 1 for all sites j becomes z^ = N5kfi- Since the system 
is translation-invariant it is more convenient to study <fi = iV _1 X)j(%) and p = 
N" 1 Ej(nj(nj - 1)). (Note that (f) = p is a constant of the motion, so this quantity 
serves principally check on our analysis.) 

Let 

Ut(U}AC})= hfD1>g'[M], (17) 



where Q' represents the exponential in Eq.(||), and let 

(A) = e- N * fvfahl>AQ'[1>A\ m , (18) 



where A is a function of the fields ip and ■0. Evaluating the derivatives in Eqs. fll5|) 
and JIB]), we have 

= iV- 1 (^ =o ) (19) 



and 



P = ^" 2 E(^-fc) • (20) 

k 
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It is convenient at this point to perform a change of variables, letting ip k = 
ip k — Ndkfi- As a result, Q' gains a factor of e Np [which cancels the prefactor in Eq. 
([Of)], and the boundary term in the argument of the exponential vanishes when we 
set z k = N5 k o • Then we have 



(A) = \Vi)Vi) AQ\^^\ 



(21) 



where 



= exp 



-N- 1 / dt'J2^-k+ / dt'd 



(22) 



and where the "interaction" is now: 



ki,k 2 ,k 3 

- 2N~ 2 u klfi i) kl %p k2 il)_ kl _ k2 . 



(23) 



Eq. (|2~2"D with Ci = defines £?o; Eq- (HH) with Q Q in place of Q defines (^4)o- 

The usual procedure at this point would be to take the continuum (small-fc) 
limit, generating a field theory for the process. While this is not our purpose in 
the present work, we note some interesting features of the model in this context. 
The first is that the interaction L contributes nothing to the quadratic part of the 
action. (This can be seen immediately from Eq.([l]): L is quartic in the fields.) 
Thus the resulting theory is nominally massless, and has no evolution at all in 
the Gaussian approximation. (Diffusion, in this model, is cooperative, requiring 
the presence of at least two particles at the same site.) The action, moreover, 
contains no parameters whatsoever; the relevant parameter p is "hidden" in the 
initial probability distribution. A further difference from continuum descriptions of 



more familiar processes such as DP p8|J19|| is that the order parameter is given by 
(ip 2 ) not (ip). 



In simulations and cluster approximations PDH231], FES clearly show a continu- 



ous phase transition between an active and an absorbing state as the parameter p 
is varied, in close analogy to more familiar examples, such as directed percolation. 
Thus at some more reduced level we might expect to find an effective field theory 
of the usual kind, with a nonzero mass, bare diffusion coefficient, and one or more 
relevant parameters. Such theories have indeed been proposed for sandpile models 
|2~6| , |2~7|1 . They include a second field (the particle density) whose evolution is cou- 
pled to that of the order parameter. We leave the systematic derivation of such a 
description, starting from the present exact action, as a topic for future work. 
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III. PERTURBATION THEORY 



A. Free expectations 



To evaluate expectations of the form of Eq. (^Tj) we write 

(A) = (AeIo dt ' c ') , 



(24) 



with the intention of expanding the exponential. Each term in the ex- 
pansion (corresponding to a diagram, as specified below), can be evalu- 
ated once we have determined the free expectations of products of fields, 
(^hin) ■ ■ ■^ l k m (T m )ip qi (K,i) ■ ■ -ijj qn (K n )} . As usual, this can be reduced to expec- 
tations of a single field and of pairs of fields. To begin, consider 



(ipk(s)) = / VipVip i> k (s) exp 



-N 



-i 



(25) 



The integrals over ipk force the condition ipk = 0, so that ipk(s) = V'fc(O) = NpSkfi. 
Evidently, the same reasoning implies that free expectation of a product of n fields 
i/j is (Np) n if all the wave vectors are zero, and is zero otherwise. (The fact that 
a single field has a nonzero expectation may appear unusual, but in fact will not 
cause any inconvenience. It is possible to make a further change of variables to 
ip' k = Tpk — Np5k t o, at the cost of introducing four new terms in £j.) 

To evaluate free expectations involving x/j, we define an operator /C 9 (t) via the 
property: 



(26) 



Integrating by parts in the exponential, we have 
6 „ 5 



5xp- k (r 



-Go 



exp 



N 



-i 



dt'^^wdt'ip-k' + pipk=o(t = 0) 



(27) 



Recalling that ijjk{t) = 0, we have 



K q {s) = -N f dr 

J s 



s 5lp- k (T) 



(28) 



Thus, 



S 



(Ms))o = ~N I V^ViP f dr \ g = 0, (29) 

where the final result is obtained via functional integration by parts. In the same 
manner we find the basic contraction (free propagator): 

(W(«)4(s))o = -N f V^ViP f dr ip k ,{u) g = N5 kl ^ k e(u-s) , (30) 

J Js 01p_ k {T) 



where denotes the step function. As is usual in this formalism, 9(0) = p5| 



The free expectation of n fields ip and n fields ip is given by the sum of all products 
of n pairwise contractions. In case there are m > n fields ip there are additional 
factors of NpS^fl associated with each uncontracted ip. 



B. Perturbative expansion 



We now have in hand all the ingredients needed to expand expectations of the 
form of Eq. fl24|). The first and second terms in £/, Eq. (p3|), correspond to vertices 
with four and three lines (a "4- vertex" and a "3- vertex", respectively). Each ipk( r ) 
must be contracted with a ip_k{ T '), where r' > r. (The required factors of ip may 
come from other vertices or from A.) We adopt a graphical notation in which fields ip 
(ip) are represented by lines entering (leaving) a vertex. All lines are oriented toward 
the left, the direction of increasing time. Fig. 1 shows the vertices associated with 
£/, and the nodes that represent and p. We refer to the latter as "sinks" since 
they have no outgoing lines. Uncontracted fields ip are called "external lines" . 

When we expand the exponential in Eq. (|24|), the n-th order term carries a factor 
1/n!, and there are n time integrations, / dt\ ■ ■ ■ J dt n , all over the interval [0, t]. We 
absorb the factor 1/n! by fixing the time-ordering t > ti > t 2 > ■ • ■ > t n > 0. This 
imposes certain restrictions on diagram topology, since a field ip must always be 
contracted with a ip at a later time. Once this ordering is imposed, the integrand 
has no further time dependence, and the time integrations yield t n /n\. 

Before formulating general rules, we study a few examples. Consider 

p = N- 2y £,(M-k<£* f£l )o- (31) 
k 

The zeroth-order term is simply: 

N- 2 Y / (^-k)o=p 2 , (32) 
k 

i.e., the initial activity for the product-Poisson distribution. At first order we have 
the diagrams (a) and (b) shown in Fig. 2. For diagram (a) there is a combinatorial 
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factor of 2, since there are two ways to contract the lines exiting the vertex with the 
fields in p; the contribution from this diagram is 

-2 [ t drN- 1 p 2 J2^k,k = ~p 2 t (33) 

In graph (b) we see that the condition that uncontracted ip fields have k = forces 
the line exiting the vertex to have k = as well. But there is then a factor of 
cjo,o = associated with this vertex, so the contribution due to graph (b) vanishes. 
In general, we can exclude diagrams in which the two lines entering a 3-vertex have 
momenta that sum to zero. 

At order t 2 we have diagrams (c) - - (f) shown in Fig. 2; (d) and (f) vanish 
for the same reason as (b). (From here on we simply ignore such diagrams.) The 
contribution of graph (c) is readily obtained: there is a combinatorial factor of 4; 
the time integrations yield t 2 /2; there is a remaining factor of N~ 2 times the square 
of the sum encountered in graph (a). The result is p 2 t 2 /2. In graph (e) there is a 
combinatorial factor of 8; its contribution is 

16{t 2 /2)p 3 N- 1 £[1 -cos 2 k\ [1 -cos k\ = ApH 2 . (34) 
k 



It is easy to see that our expansion conserves the particle density. With a one- 
line sink in place of the two-line sink corresponding to p, there will never be enough 
factors of ip to contract with all the factors, if use only four-line vertices. On the 
other hand, if we contract the sink with a three-line vertex, there will be a factor of 
Wo,o = 0. Thus <p{t) = 0(0) =p. 



C. Diagram rules 

In the analysis that follows it will be convenient to employ the Laplace transform; 
the factor t n /n\ then becomes l/s n+1 , where s is the transform variable. Each 
diagram in the series for p(s) carries a factor of tv^ -3 ™ 4-2 ™ 3-2 ; where L is the number 
of lines (including external lines), and n 3 and n 4 are the numbers of 3- vertices 
and 4- vertices, respectively. There is exactly one factor of iV -1 for each free wave 
vector sum; from here on, we simply associate such a factor with each sum. We 
can formalize the foregoing discussion into a set of rules for finding the ra-th order 
contribution to the order parameter p(s): 

1) Draw all connected diagrams with n vertices, and a two-line sink to the left 
of all vertices. It is permissible for a line entering a node to be uncontracted (each 
external line carries a factor p and must have momentum zero), but each line exiting 
a vertex (J) must be contracted with a vertex (i < j) to the left; there is a factor 
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of 8k' -k associated with each internal line. (Here k is the wave vector of the line 
exiting vertex j and k' the wave vector entering vertex i.) Given the restriction 
noted above, that the sum of the momenta entering a 3-vertex cannot be zero, the 
rightmost vertex of any diagram must be a 4-vertex. We refer to this (and any 
other) 4-vertex with two external lines as a source point. 

2) For each graph there is an overall factor of 1/s, and a combinatorial factor count- 
ing the number of ways of realizing the contractions. In the series for p, the combi- 
natorial factor is the product of a factor Cy, associated with the choice of lines at 
each vertex (for a fixed set of connections between vertices), and a factor Cl giving 
the number of choices of, and connections between vertices, consistent with a given 
diagram. It is easy to show that Cy = 2 Q , where Q = 1 + n 3 + 2n 4 — £ — f, with 
I the number of simple loops, and / the number of source points. By a simple loop 
we mean a pair of vertices directly connected by two lines, as in Fig. 2 (a). The 
connection factor Cl is unity for diagrams of three or fewer vertices, but can take 
nontrivial values for n > 4. Examples are discussed below. 

3) Each 3-vertex carries a factor of — 2s~ 1 Uk t0 = — 2s~ 1 [l — cos A;], where k is the 
momentum leaving the vertex. Each 4-vertex carries a factor of — s -1 ^^, where 
k\ and ki are the momenta of the lines exiting the vertex. 

4) After taking into account all of the 5-functions associated with propagators, the 
remaining free wave vectors are summed over the first Brillouin zone, Eq. (^). 

We close this subsection with a discussion of the combinatorial factor Cl- As 
noted, Cl represents the number of distinct choices of vertices, and of connections 
between vertices, consistent with a given diagram topology. Evidently, diagrams 
with n vertices arise when we expand the product Cj^ by "choice of ver- 

tices" we mean whether the 3-vertex, or the 4-vertex, associated with each Li^ is 
selected. The n-th vertex is, as noted, always a 4-vertex. Of course, if the remaining 
n—1 vertices of a diagram are all of the same kind, then only one choice exists. 
In most cases, exchanging the positions a 3-vertex and a 4-vertex yields a different 
diagram, but this is not always so. Consider, for example, diagram (a) shown in 
Fig. 3. Let us refer to the i-th vertex in a given diagram (in order of decreasing 
time) as V{, with Vo denoting the sink. In this diagram V\ must be a 3-vertex, and 
V4 a 4-vertex, but we are free to choose between V2 and V3 as the other 3-vertex. 
Thus this diagram appears twice in the expansion of the activity, so that Cl = 2 in 
this case. (That Cl = 2, and not more, rests on the fact that, given the choice of 
vertices, there is only one way of connecting them to yield the desired topology.) 

Next we consider the possibility of different connections among a fixed set of 
vertices. We use to denote a link between Vi and Vj. Consider diagrams (b) 
and (c) of Fig. 3. For diagram (b), the set of connections between vertices is unique: 
[0,1], [0,4], [1,2], [2,3], and [3,4], so Cl = 1 for this diagram. For diagram (c), by 
contrast, there are several different sets of connections that yield the same topology. 
Evidently, the links between V\ and V , and between V4 and V3, are obrigatory. But 
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this leaves open the question of which vertex the other outgoing line of V4 links to, 
and of which vertex is linked to the other line entering the sink. One readily verifies 
that the following possibilities exist: 



i) [0,1], [0,2], [1,4], [2,3], [3,4]; 

ii) [0,1], [0,2], [1,3], [2,4], [3,4]; 

iii) [0,1], [0,3], [1,2], [2,4], [3,4]. 

Thus Cl = 3 for diagram (c). Consider, finally, diagram (d) of Fig. 3. If V\ and V 2 
are both 3-vertices, there are two possible sets of connections: 

i) [0,1], [0,2], [1,3] 2 , [2,4], [3,4]; 

ii) [0,1], [0,2], [2,3] 2 , [1,4], [3,4]. 

In addition, this diagram can be realized (with a unique set of connections), with 
V\ and V 3 as the 3-vertices, so that (7l = 3 in this case. 



D. Examples and some general results 



There are seven diagrams at third order (Fig. 4). Diagram (a), which continues 
the series whose first two terms are (a) and (c) of Fig. 2, is readily shown to yield 
—p 2 /s 4 . Applying the rules to diagram (b) in Fig. 4, we find: 

16(-2)iV-V^^[l-cos 2 A;]^[l-cos 2 A; , ][l-cosA; / ] = -^p 3 . (35) 
s k k> s 

The contribution of diagram (c) is identical. For diagram (d) we have 

16(-4)A^y^]T[l-cos 2 A;][l-cosA;] 2 = -—^ . (36) 

s k s 



Diagram (e) makes the same contribution. Similar analysis yields —18p 3 /s 4 for 
diagram (f) and — 32p 3 /s 4 for (g). Collecting results, we have 



2 + 3 



p=p 2_ p 2 t+ ii_^ + 8p) _ii_ 



l + 66p + 80p 2 +0(t 4 ) 



(37) 



In evaluating wave vector sums the following observations are helpful. First, note 
that any wave vector sum iV -1 J2k can be written as (27r) _1 JZ n dk. The analysis is 
facilitated by use of the identity: 



12 



In = ^E(l"C°s 2 fc) [l-oos fcf = 4 ■ 2"|^|j ■ (38) 



[This can be shown by writing I n as an integral, and letting u = l — cos A; so that 

1 r 



/ n = _ (* duu n+1/2 V2^i . 

IT JO 



Integrating by parts one readily finds that 

/ -2 2n+1 I 
n ~ 2n+4 n - 1 ' 

leading to Eq (|38"D , since Jo = 1/2. For convenience we note: l\ = 1/2; I 2 = 5/8; 
J3 = 7/8; and J4 = 21/16.] Useful identitites in the ^-dimensional case are: 



and 



]TA d (k)A d (q-k) = ^A d ( q ). 
k - 



Two general properties of the expansion are readily verified. First, the leading 
term at each order is p 2 (—t) n /n\. This set of contributions sums to p 2 e~ f , describing 
the relaxation of the initial density of active sites. In d dimensions this becomes 
p 2 exp[— (2d— l)t/d]. The decay rate, (2d—l)/d, represents twice the probability that 
the two particles released by a toppling site move to distinct neighbors, the factor 
of two representing the toppling rate for a site with exactly two particles. 

A second observation is that the highest power of p appearing at order t n is p n+1 , 
coming from diagrams with one source point and n — 1 external lines terminating 
at 3-vertices. The coefficient of this term can be found as follows. Note that the 
relevant diagrams consist of (n—1) 3-vertices attached to the lines of the simplest 
graph (Fig. 2a). The number of realizations of such topologies (that is, the sum of 
the Cl over all diagrams of this kind), is given by one half the number of distinct 
sequences of r symbols A and n — r — 1 symbols B (A and B corresponding to 3- 
vertices inserted in one or the other line), summed over r from r = to n — 1. (The 
factor of one half is required because a pair of sequences that merely exchange A's 
and .B's in fact represent the same diagram.) The number of such realizations is: 
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- VV n ~ l \ = 2 n ~ 2 
2 h V ' J 

Each diagram of this kind has a combinatorial factor Cy = 2 n+l , and an additional 
factor of (— l) n 2 n ~ 1 / ra _ 1 coming from the 3- vertices and the wave vector sum. Using 
our result for I n , we find that this contribution is: 

, n ng4n-i (2n-l)!!p w+1 t" 
1 1 (2n + 2)!! n\ ' 

for n > 2. Call the sum of such contributions p ma x(t); it can be evaluated in closed 
form by noting that 



pi /■*,, 2f ~ (-8p*) n (l-cosfc) 



Pmax(t) = 7— / dk(l-COS 2 k) V 
4 Z7T J-7T n 



n-1 



n=2 nl 



= V - {e- 8 ^ [I (8pt) + h(8pt)} - 1} + p 2 t, (39) 

where I v denotes the modified Bessel function. The factors of e~ 8pt I u (8pt) imply 
slow relaxation (~ typical of diffusive processes. 



IV. DIAGRAMMATIC ANALYSIS 



As is evident from the foregoing discussion, many diagrams represent simple 
variations of those appearing at some lower order. The calculation will be simplfied 
if we can identify classes of related diagrams. To begin we enumerate some ways in 
which diagrams (g') with v + 1 vertices can be formed from a f-vertex diagram g. 

a) Insert a 3-vertex into any internal line; g' has one more external line than does 
g. An example is the generation of diagram (e), starting from (a) in Fig. 2. 

b) Replace any 4-vertex by a "4-loop" , consisting of a pair of 4-vertices joined by 
two lines. Diagram (c) (Fig. 2) is generated from (a) by this process. 

c) Replace any 3-vertex V by a 4-vertex, contracting the new outgoing line with 
a new 3-vertex inserted in an existing internal line. In case the contraction is to 
a 3-vertex V on the same line as V, and there are no other vertices on this line, 
between V and V, we call the resulting structure a "3-loop". Starting with diagram 
(e) in Fig. 2, we may generate (f) and (g) (Fig. 4) by this procedure. Diagram (f) 
has a 3-loop but (g) does not, since in this case V does not lie on the same line as 
as V. 

d) Insert a 4-vertex into an internal line, contracting one of its outgoing lines to an 
existing external line. (The vertex associated with the latter must lie to the left of 
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the new 4-vertex.) Graphs (f) and (g) (Fig. 4) are formed from (e) (Fig. 2) by this 
means. 

e) Add a new 4-vertex to g by joining its outgoing lines to two external lines in g. 
Consider again (e) in Fig. 2; if we join two external lines (one of them associated 
with the 3-vertex) to a new 4-vertex, the result (after some redrawing) is (g) of Fig. 
4. 

Now consider the inverse of operations a) - c). Given a graph g with two or more 
vertices: 

a') Remove all 3-vertices bearing an external line. 

b') "Collapse" any 4- loop to a single 4-vertex. 

c') Collapse any 3-loop to a 3-vertex. 

The systematic application of these steps, until no further subtractions are pos- 
sible, will be called reduction of a diagram. Let us define an articulation point as 
a 4-vertex (other than a source point), which, if removed, will leave the diagram 
disconnected. (Such a vertex has, by necessity, total momentum zero entering it.) 
In Fig. 2 only graph (c) has an articulation point, while graphs (a), (b) and (c) of 
Fig. 4 all possess one or more such point. 

We may now define an irreducible diagram (IRD) as one free of 3-vertices bearing 
external lines, free of collapsible loops, and having no articulation points. The 
simplest irreducible diagram is (a) of Fig. 2; we shall refer to it as the one-loop 
IRD. At each order only a small minority of the diagrams are irreducible. At second 
order there are no IRDs; the only one at third order is (g) in Fig. 4. The IRDs with 
four vertices are shown in Fig. 5; there are 23 IRDs with five vertices. 

In addition to the local modifications (a - e) described above, we can also form 
new diagrams via composition. Given two diagrams g and g' we form the composite 
diagram gg' by: 

i) Deleting the two lines entering one of the source points r of g; 

ii) Identifying the sink in g' with r in g. 

(If g and g' are distinct we can also form g'g.) Graph (b) of Fig. 4 is composed, in 
this sense, from (a) and (e) of Fig. 2. Composite graphs are reducible by definition 
(the two components are joined at an articulation point.) When we have removed 
all 3-vertices with external lines, and collapsed all collapsible loops in a diagram, 
the result is either an IRD or a composite of such diagrams. 

The expansion up to order t v may be organized as follows. First we identify all of 
the IRDs with v or fewer vertices. For each such IRD, we evaluate its contribution 
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and that of all diagrams with < v vertices that are directly reducible to it. There 
remain the composite diagrams. Let C g (s) be the contribution to p(s) due to diagram 
g. From the definition of the composition process, we see that the contribution of a 
composite diagram gg' to the activity is: 



gg 



P 



2 CgCgl. 



(40) 



We make use of this relation as follows. For each IRD g*, define 



F g *(s) 



g'yg* 



(41) 



as s/p 2 times the sum of contributions due to g* and all graphs reducible to it (we 
use g' y g* to mean u g' is reducible to g* n ). Now define 



g* 



(42) 



and 



(43) 



where the sum is over all IRDs, and and in the second sum n g * is the number of 
source points in g*. Each source point serves as a possible attachment site for the 
preceding element. Then the sum of (i) all irreducible diagrams, (ii) all diagrams 
reducible to an IRD, and (iii) all linear composite diagrams is 



Pl{s) 



F(s) 



s 1 - G(s) 



(44) 



Missing from pi are branching composite diagrams. These are composite diagrams 
in which at least one component has two or more diagrams attached to its source 
points. Since diagrams with more than one source point have at least four vertices, 
branching diagrams must have at least six. Thus p(s) = pi{s) to 0(s~ & ). Evaluating 
F to and G to this order, we obtain 

p = p -p t + ——(1+Sp) - 



6 



.2+4 



+ 



p z t 



1 + 66p + 8(V 
1 + 442p + 2076p 2 + 896p 3 

1 + 2842p + 35396p 2 + 52240p 3 + 10752p 4 l + 0{f ) . 



(45) 
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In d dimensions the result is: 



p 2 t 2 




p 2 t 3 
6 


7 3 + ^ (80rf 4 -56rf 3 + 12rf 2 -6d+3) + ^jf- (8d 3 -6d+3) 




p 2 t 4 
+ 4! 


7 4 + ^ (ll52d 5 -960d 4 + 336d 3 -128d 2 + 36d+6) 




d 5 v 


960d 5 368d 4 112d 3 48d 2 + 108d 2l) + J (8d 2 (d+l) 9(2d 1)) 


+ 0(^ 5 ) , 


(46) 



where 7 = (2cf— l)/d. 



V. COMPARISON WITH SIMULATION 



Although the series we have derived are too short to yield reliable predictions for 
the activity p(t) at long times, it is of interest to consider a preliminary comparison 
with simulation results, as a check on our analysis. We have simulated the one- 
dimensional model on lattices of L =1000 and 2000 sites, with periodic boundaries; 
initially N=pL particles are placed randomly on the lattice, generating a product- 
Poisson distribution of occupation numbers rij. We average over 5 x 10 5 independent 
realizations of the process. 

A typical result is shown in Fig. 6, for p= 1/2. Since this is well below the critical 
value, p(t) — > as t — > 00. (Note that the average is over all realizations, including 
those that have fallen into the absorbing state.) The simulation result (data points 
merging to the bold line) decays quite rapidly at first, and then approaches zero 
more slowly; a systematic study ||T| reveals that the approach to p = is best 



characterized as a stretched exponential, p ~ exp[— at?] with {3 ~ 1/2. 

Two theoretical curves derived using the series, Eq. ([45|) , are shown in Fig. 
6. The upper one is generated by transforming the time series to the variable 
y = [1 — e~ bt ]/b, and then constructing the [3,2] Pade approximant to the y-series. 
For b in the range 0.25 to 0.4 we find good agreement at short times; for b = 0.35 
(shown here) there is excellent agreement for t < 10, but the theoretical prediction 
attains a constant, nonzero value thereafter. 



Since pit) must decay to zero in this case, it is interesting to study the series 
for dlnip/p 2 )/dt. The latter remains negative for large t, so that the activity itself 
decays to zero. The lower curve in Fig. 6 is obtained by transforming the derivative- 
log series to the variable z — 1 — (l + fot) -1 / 3 (with 6 = 2.7), and forming the [2,2] 
approximant to the ^-series, which is then used to evaluate ln(p/p 2 ) via numerical 
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integration. The activity does indeed decay to zero, but too rapidly. We obtain 
very similar results using b in the range of 2 to 3, and using other transformations, 
for example y given above, or a modified z with the exponent 1/2 in place of 1/3. 
After examining various further transformations and analyses, we conclude that the 
present series is simply too short to yield useful predictions for long times. 

It should be noted that the apparent radius of convergence of the original series 
(estimated from ratios of successive coefficients), is about t < 1/4 for p = 1/2; ev- 
idently, the transformation and approximant analysis greatly extends the utility of 
the series. (For p = l, the radius of convergence is about t < 1/6; the transformed 
series agrees well with simulations for t > 2.) Thus, despite the divergence between 
simulation and theory for times greater than about 10, we regard the present com- 
parison, based as it is on a very short series, as confirming the validity of our analysis. 
Since the approach to the stationary state is neither exponential nor power-law, we 
anticipate that a fair number of terms will be required to represent it. Once ex- 
tended series are available, transformations that are consistent with the asymptotic 
behavior should yield good numerical estimates. 



VI. DISCUSSION 



We have derived an exact path-integral representation for the dynamics of a 
stochastic sandpile model similar to Manna's sandpile, in a version with strictly 
conserved particle number. (The key differences are: (1) We analyze a continuous- 
time process, whereas Manna's model is defined in discrete time; and, (2) The 
toppling rate for a site with n particles is n(n— 1) rather than uniform, for all sites 
with n > 2.) We have argued that such minor differences should not affect scaling 
properties. The path-integral mapping yields an effective action that is nominally 
massless and parameter- free, the relevant parameter p (the particle density), appear- 
ing in the initial distribution not the evolution operator. This conclusion applies 
generally to non-dissipative sandpiles. 



In contrast to phenomenological field theories of sandpiles p6 , P?[ , in which the 
order parameter density is coupled to a second slowly relaxing variable (the particle 
density), the exact action derived in the present work involves but a single field, 0, 
whose expectation is the particle density. The order parameter is given by (ip 2 )- 
The derivation of a phenomenological description along the lines of Refs. |26| and 



27] , starting from the exact action, remains an open question. 



We develop a diagrammatic perturbation theory, leading to an expansion for the 
activity density in powers of time, where the coefficients are polynomials in p. We 
point out certain general properties of the expansion: its limiting forms for very 
small, and very large p, and the slow, diffusion- like nature of the relaxation. The 
series provides an important check on other theoretical methods, and on simulation 
results. While it appears to have a limited radius of convergence, experience with 
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similar series, for the contact process |3T] and for random sequential adsorption 



42|,32], indicates that under suitable transformation of variables and Pade approx- 
imant analysis, useful predictions can be obtained. We present some preliminary 
numerical evidence for this assertion. Thus, if the series can be extended, it should 
be possible to investigate scaling properties of the sandpile. We intend to pursue 
such extensions and analyses in the near future. 
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FIGURE CAPTIONS 



Fig. 1. Vertices and sinks in Ci, p, and 0. 

Fig. 2. One- and two-vertex diagrams in the expansion of the order parameter. 
Fig. 3. Diagrams analyzed in the discussion in the combinatorial factor Cl- 
Fig. 4. Diagrams with three vertices. Only (g) is irreducible. 
Fig. 5. Irreducible diagrams with four vertices. 

Fig. 6. Activity in the sandpile with p= 1/2. The points (merging to a bold line) are 
simulation results; the curves above and below are series-based predictions discussed 
in the text. 
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